# SimesJahn.R
# Main Results
# "The Consequences of Medicaid Expansion Under the Affordable Care Act for Police Arrests"
# Last updated by JS: 12/9/21

library(dplyr)
library(ggplot2)
library(sandwich)
library(cowplot)

load("SimesJahn.RData")

##-----------------------------------------------------------------------
## Descriptive Statistics Table
##-----------------------------------------------------------------------

# Table 2
tab2 <- analysis[analysis$Year<2014,]
tab2 %>% summarise(all=mean(arr,na.rm=TRUE),violent=mean(vlr,na.rm=TRUE),drug=mean(drr,na.rm=TRUE), 
                   lowlevel=mean(llr,na.rm=TRUE),black=mean(black,na.rm=TRUE),mage=mean(medage,na.rm=TRUE),
                   pov=mean(poverty,na.rm=TRUE),unemp=mean(unemp,na.rm=TRUE),urb=mean(urban,na.rm=TRUE),
                   welf=mean(exsocial,na.rm=TRUE),edu=mean(exeduc,na.rm=TRUE),opioidd=mean(death, na.rm=TRUE))

tab2 %>% group_by(treatment) %>% summarise(all=mean(arr,na.rm=TRUE),violent=mean(vlr,na.rm=TRUE),drug=mean(drr,na.rm=TRUE), 
                                           lowlevel=mean(llr,na.rm=TRUE),black=mean(black,na.rm=TRUE),mage=mean(medage,na.rm=TRUE),
                                           pov=mean(poverty,na.rm=TRUE),unemp=mean(unemp,na.rm=TRUE),urb=mean(urban,na.rm=TRUE),
                                           welf=mean(exsocial,na.rm=TRUE),edu=mean(exeduc,na.rm=TRUE),opioidd=mean(death, na.rm=TRUE))

##-------------------------------------------------------------------------
## Models and Sandwich SEs
##-------------------------------------------------------------------------
library(MASS)
m1<-glm.nb(all ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
           + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath
           + offset(log(countypop)), data=analysis)

m2<-glm.nb(violent ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
           + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath
           + offset(log(countypop)), data=analysis)

m3<-glm.nb(drug ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
           + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath
           + offset(log(countypop)), data=analysis)

m4<-glm.nb(lowlevel ~ as.factor(state) + as.factor(Year)+ as.factor(expansion)
           + medage + black + unemp + poverty + urbanicity + exsocial + exeduc + ldeath
           + offset(log(countypop)), data=analysis)

ses<- function(model) { 
  cov <- vcovHC(model, type = "HC0")
  std.err <- sqrt(diag(cov))
  q.val <- qnorm(0.975)
  r.est <- as.data.frame(cbind(
    Estimate = coef(model)
    , "Robust SE" = std.err
    , z = (coef(model)/std.err)
    , "Pr(>|z|) "= 2 * pnorm(abs(coef(model)/std.err), lower.tail = FALSE)
    , LL = coef(model) - q.val  * std.err
    , UL = coef(model) + q.val  * std.err
  ))
  return(r.est)
}
mod1<-ses(m1)[52:54,]
mod2<-ses(m2)[52:54,]
mod3<-ses(m3)[52:54,]
mod4<-ses(m4)[52:54,]

mods<-rbind(mod1,mod2,mod3,mod4)
mods$Estimate<-exp(mods$Estimate)
mods$LL<-exp(mods$LL)
mods$UL<-exp(mods$UL)
mods$model<-c(rep("All",3),rep("Violent",3),rep("Drug",3),rep("Low Level",3))
mods$year<-rep(c(1,2,3),4)

##-------------------------------------------------------------------------
## Fig 1 and Fig 2
##-------------------------------------------------------------------------

# Fig 1: Time series plot of all arrests and drug arrests

plot<-analysis[analysis$Expansion15!=1 & analysis$Expansion16!=1,]

plotdata<-plot%>% 
  group_by(Year,Expansion14) %>%
  summarise(drug=mean(drr, na.rm=TRUE), arrests=mean(arr, na.rm=TRUE))
plotdata<-plotdata[!(is.na(plotdata$arrests)==TRUE),]
plotdata$Expansion14<-ifelse(plotdata$Expansion14==1, "Expansion States", "Non-Expansion States")

#pdf("fig1.pdf")
pa<-ggplot(plotdata, aes(x=Year, y=arrests, color=Expansion14))+
  geom_line(aes(linetype=Expansion14))+geom_point(aes(shape=Expansion14, fill=Expansion14),size=3)+
  scale_linetype_manual(values=c("dotdash", "solid"))+
  scale_shape_manual(values=c(1, 16))+
  geom_vline(xintercept=2014,linetype="dashed", color="grey50")+
  ylab("Arrests per 100,000 Population")+ xlab(" ")+ ggtitle("All")+
  theme(legend.position="bottom", legend.title = element_blank(),legend.direction = "horizontal",legend.key = element_rect(colour = "transparent", fill = "white"),
        plot.title = element_text(hjust = 0.5),text = element_text(size=15,colour = "black"),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))+
  scale_color_manual(values=c("#FC4E07","#00AFBB"))
pb<-ggplot(plotdata, aes(x=Year, y=drug, color=Expansion14))+
  geom_line(aes(linetype=Expansion14))+geom_point(aes(shape=Expansion14, fill=Expansion14),size=3)+
  scale_linetype_manual(values=c("dotdash", "solid"))+
  scale_shape_manual(values=c(1, 16))+
  geom_vline(xintercept=2014,linetype="dashed", color="grey50")+
  ylab("Arrests per 100,000 Population")+ xlab(" ")+ ggtitle("Drug")+
  theme(legend.position="bottom", legend.title = element_blank(),legend.direction = "horizontal",legend.key = element_rect(colour = "transparent", fill = "white"),
        plot.title = element_text(hjust = 0.5),text = element_text(size=15,colour = "black"),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))+
  scale_color_manual(values=c("#FC4E07","#00AFBB"))
p<-plot_grid(pa+ theme(legend.position="none"),
             pb, nrow=2,  align="v")
print(p)
#dev.off()

# Figure 2: DID Results

mods$Change<-(mods$Estimate-1)*100
mods$CUL<-(mods$UL-1)*100
mods$CLL<-(mods$LL-1)*100

#pdf("fig2.pdf")
p1<-ggplot(mods[mods$model=="All",], aes(x=year, y=Change))+
  geom_point(color="#00AFBB",size=3)+
  geom_line(color="#00AFBB")+
  geom_errorbar(aes(ymin=CLL, ymax=CUL),color="#00AFBB", width=0.2)+
  coord_cartesian(ylim=c(-51,4))+ 
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ylab("% Change")+ggtitle("All")+
  scale_x_discrete(limits=c("Year 1", "Year 2", "Year 3"))+
  theme(legend.position="bottom",legend.title = element_blank(),legend.direction = "vertical",legend.key = element_rect(colour = "transparent", fill = "white"),
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black"),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), axis.title.x=element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p2<-ggplot(mods[mods$model=="Violent",], aes(x=year, y=Change))+
  geom_point(color="#00AFBB",size=3)+
  geom_line(color="#00AFBB")+
  geom_errorbar(aes(ymin=CLL, ymax=CUL),color="#00AFBB", width=0.2)+
  coord_cartesian(ylim=c(-51,4))+ 
  #scale_color_manual(values="#00AFBB")+
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ylab("% Change")+ggtitle("Violent")+
  scale_x_discrete(limits=c("Year 1", "Year 2", "Year 3"))+
  theme(legend.position="bottom",legend.title = element_blank(),legend.direction = "vertical",legend.key = element_rect(colour = "transparent", fill = "white"),
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black"),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), axis.title.x=element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p3<-ggplot(mods[mods$model=="Drug",], aes(x=year, y=Change))+
  geom_point(color="#00AFBB",size=3)+
  geom_line(color="#00AFBB")+
  geom_errorbar(aes(ymin=CLL, ymax=CUL),color="#00AFBB", width=0.2)+
  coord_cartesian(ylim=c(-51,4))+ 
  #scale_color_manual(values="#00AFBB")+
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ylab("% Change")+ggtitle("Drug")+
  scale_x_discrete(limits=c("Year 1", "Year 2", "Year 3"))+
  theme(legend.position="bottom",legend.title = element_blank(),legend.direction = "vertical",legend.key = element_rect(colour = "transparent", fill = "white"),
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black"),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), axis.title.x=element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p4<-ggplot(mods[mods$model=="Low Level",], aes(x=year, y=Change))+
  geom_point(color="#00AFBB",size=3)+
  geom_line(color="#00AFBB")+
  geom_errorbar(aes(ymin=CLL, ymax=CUL),color="#00AFBB", width=0.2)+
  coord_cartesian(ylim=c(-51,4))+ 
  #scale_color_manual(values="#00AFBB")+
  geom_hline(yintercept=0, linetype="dashed", color="grey50")+
  ylab("% Change")+ggtitle("Low Level")+
  scale_x_discrete(limits=c("Year 1", "Year 2", "Year 3"))+
  theme(legend.position="bottom",legend.title = element_blank(),legend.direction = "vertical",legend.key = element_rect(colour = "transparent", fill = "white"),
        plot.title = element_text(hjust = 0.5),text = element_text(colour = "black"),
        panel.grid.minor = element_blank(),panel.grid.major = element_blank(), axis.title.x=element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"))
p<-plot_grid(p1+ theme(legend.position="none"),
             p2+ theme(legend.position="none"),
             p3+ theme(legend.position="none"), 
             p4+ theme(legend.position="none"), nrow=2,  align="h")
print(p)
#dev.off()